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Abstract. By using the Hugoniot curve in detonics as a Riemann invariant of a velocity-pressure 
model, we get a conservative hyperbolic system similar to the Euler equations. The only differences 
are the larger value of the adiabatic constant (7 — 8.678 instead of 1.4 for gas dynamics) and the 
mass density replaced by a strain density depending on the pressure. The model is not homogeneous 
since it involves a gravity and a friction term. After the seismic wave reaches up the bottom of 
the ocean, one gets a pressure wave propagating toward the surface, which is made of a frontal 
shock wave followed by a regular decreasing profile. Since this regular profile propagates faster than 
the frontal shock waves, the amplitude of the pressure wave is strongly reduced when reaching the 
surface. Only in the case of a strong earth tremor the residual pressure wave is still sufficient to 
generate a water elevation with a sufficient wavelengths enable to propagate as a Saint Venant water 
wave and to become a tsunami when reaching the shore. We describe the construction of the model 
and the computation of the wave profile and discuss about the formation or not of a wave. 

We propose a model using a constant mass density because variable mass density models are often 
unstable, since a tiny variation of the density allways causes a large variation of the pressure. In 
the model presented here, the transported variable is not the mass density but a new variable, 
called the strain density, which has the same properties of conservation without this drawback of 
numerical unstability. 

The first section deals with the construction of the velocity-pressure model and the new con- 
servative variable of strain density is designed in the second section. The new non homogeneous 
Euler-like model is then studied in the third section and one dimension numerical computations 
of the profile of the wave are reported in the fourth section, for different values of the friction 
coefficient. We conclude by some discussions about the emergence of a tsunami wave or not. 

1. The velocity-pressure model 

The Hugoniot curves (see [6]) correspond to the linkage between the velocity w and the pressure p 
in a shock wave travelling through a material after an impact. By starting from a position at rest, 
the algebraic equation of such curves has the form 

p = po { Co w + So w'^) , 

where po is the mass density of the material, cq is its speed of sound, and 5*0 is a dimensionless 
constant retated to this material, obtained from experiments. For water, po = 1000 kg m~^, cq = 
1647 m and So = 1.921. By solving with respect to w , we get 



w{p) = ± ao ( a/1 + Pop - 1 
and Po = ^ = 2.84 10"^ whoe 

Po ^0 

We notice the relation ao Po po Cq = 2 . The values of these parameters come from [5]. 



where ao = 7t%- = 429 m s ^ and do = = 2.84 10 ^ whose unit is the inverse of an energy. 
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A general one dimension velocity-pressure model, with a constant mass density po; is made of 
a dynamics equation of the form 

(1.1) po (wt + ww^) + = 

where t is the time and z is the position along a vertical upwards oriented axis (with 2; = at sea 
level) , and a Hooke law of the form 

(1.2) Pt + w Pz + Po cipf = 

where c{p) > stands for the pressure depending wave velocity, to be identified. These two equations 

compose a hyperbolic system whose Riemann invariants have the form w'(p) — ± 1-— . By 

Po c{p) 

derivating the expression of w{p) from the Hugoniot curves, one gets 



w'(p) = ± 



Identifying the two expressions and using ao/^oPoCo = 2 lead to the formula 

c{p) = Co V 1 + /3oP ■ 

This reads like a state law for our velocity pressure model. 

Up to now we were only concerned with homogeneous equations since Riemann invariants only 
exists in this case. We also have to take in account the gravity effects, since the pressure is increasing 
with the depth of water, and some friction effect since stillness is a stable configuration. We use 
a friction term of the Strickler type as usual in hydraulics, but other choices are possible and will 
lead to similar results. That way, the dynamics equation is replaced by 

(1.3) Po (wt + wwz) + Pz + Po g + k \w\w = , 

where g is the gravity constant and k the friction parameter. The Hooke law and the state law 
are unchanged, and we take g = 9.8 m s~^m the numerical experiments. The size of the friction 
parameter is a priori unknown and will be discussed later. 

2. The conservative strain density 
We look for a quantity q — q{p) satisfying the transport equation 
(2.1) qt + {qw)^ = . 

Since q depends on p we get 

q'{p) {Pt + w Pz) + qip) = , 
to be compared with the Hooke law. We get 

q\p) ^ 1 ^ 1 
q{p) Po c(p)2 Po cl (l + /3op) ' 
which is easily solved and gives 



qip) = (1+ l3op)^°''°'^ = (1 + /3op)^=o , 
for the choice g(0) ~ 1 . The wave velocity can be written as a function of the strain density q as 

c{q) = Co . 
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We check now the conservation of the momentum m = qw . We get 

q k q 

rrit + [rnw)^ + — Pz + gg + \w\w = , 

Po Po 

where 

Q 2 Co 2^ 2 2^ , .2 

— Pz = 3 q °'o Qz = Cq q °o q^ = c[q) q^ • 

Po "0 Po Po 

We introduce a pressure term, standing as a strain pressure, 

2^+1 2 

= cl^^ = ^q^\ njith^o = 2^ + l{= 8.678) 
and we get the conservative equation for the momentum 

k 

(2.2) rrit + {mw + P{q)).^ + + — q\w\w = 0. 

Po 

The system made of (12.1 p (12.21) has the form of the well known Euler equations for gas dynamics 
with a larger adiabatic coefficient 70 = 8.678 instead of the usual value 1.4 for gases, and can be 
handled in the same way, especially for the shock waves. The same Rankine Hugoniot condition is 
valid, connecting the velocity of a shock wave to the two states {qi,wi) and {q2,W2) by 



(^2 3) ^/(^^) = W1+W2 ^ qi + q2 / 1 ^^(^2) - P{qi] 



qiq2 q2 - q\ 



CQ 

Recalling that g "o = 1 + /?Qp, we can use the atmospheric pressue Pa at the surface [z = 0) and get 



3. The profile of the strain wave 
We first compute the state at rest. In case of stillness the equations (12. ip and (12. 2p reduce to 

qt = , Co g^^"^ qz + 9 = , 
since qy^O. By denoting q = qo{z) the strain density at rest, the integration gives 

°^ ° go(-2)^"° + 9 z = Constant . 
; that g "0 = 1 + /?Qp, we cai 

"'^ ^° ^° (p-Pa) + 5-2; = , that is p = pa-po9z 01 qo{z) = ( 1 + /3o (pa- po9z))^ 

which corresponds to the geostrophic equilibrium state. 

To compute the strain density profile we use the deviation variable f] = q — q^ and look for 
linkage of the form m = A rj — i? , as in any q — m— system with a source term (see [4], or 
annex below in Section 5), where A and B are constant. Since stillness is reached for rj = 0, we 
have 5 = 0. Besides, since go does not depend on t, we have 

r]t + At]^ = , 

which means that A corresponds to the wave velocity, that we name the reference velocity, corre- 
sponding to a reference state (gre/, qrefWref) such that A = Wref + c{qref)- We compute, for q > qo 
which is always expected, 

^ A q~ qo , ,2 A 9qo 

w = — = A , nit = Arjt = -A r]^ , = Ar]^ , q^ = Vz 



ciqof 
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to be introduce into (12.21) which becomes 



A 



~ % Y 
1 



Vz +gq- 



This equation reduces to 

Vz + gq + 



a4 

r. 



gqo 
c(go)' 



A' 



c(go)^ 



2 / g -go 
q 



A 



2iq 



go. 



+ — T (g - q^y 

po q 



. 



+ ^^^{q-qoY 
Po q 



. 



By multiplying by and using c = cq q°'o we get 



(3.1) (g^c^-ggM^) ry, + 



v 



x2^-l\ 

qo) 

J 



+ 



+ [q-qo) 



c{qo 



Po 



Since q = rj + qq, this is a differential equation which can be integrated by using standard numerical 
methods. An increasing profile is expected. Since the reference velocity A is far larger than c{q) the 
coefficient g^c^ — gg^^ of rj^ is always negative in practice. The two last terms are always positive, 



and the friction term is the predominant one. The term gq^ 



1 - 

QO 



2^-1 
q \ "0 



is always negative 



and is always balanced by the friction term when the friction coefficent k is not too small. 

The value of A is determined by the strength of the sismic wave at the bottom of the ocean, 
whose depth is denoted Zf. This corresponds to a reference state {qref,^ref)- We have 



A 



Wref + Cref , rUref = qrefWref = A {qref - qo{z f)) , C^e/ 



£0. 

Co qref 



We get 



Qref 



A 



and w. 



A 1- 



-re/ 



qojzf) 

qref 

Qref 



-^ref ) 



Whichs gives A = Cref ^7-^ ana Wref — ^rer \ I \ 

Now we can compute the profile of the stain wave as the solution of (13. ip . 
Figure 1 presents a series of numerical computation tests using the reference value g^e/ = 1.1296, 
for a depth Zf = 3700 meters, that is an increasing of about 1500% above the natural pressure on 
the bottom of the ocean. 



The velocity field Wref increases from a few meters per second near the bottom to more than 
300 meters per second near the front shock wave drawn here, which is here an hypothetic one. The 
tests performed with too small friction coefficients (observed here for k < 0.15) yield decreasing 
profiles, as expected from the remark above about the size of the friction term. The real front shock 
wave will progress more slowly than the strain wave, and its amplitude will decrease rapidly, and 
the effective values of the velocity field Wref will be strongly reduced near the front schock wave. 
The front shock wave connects the geostrophic equilibrium state go, with the velocity w = 0, to a 

value g on the strain wave, with the velocity w = A (^1 — ^ j • The Rankine Hugoniot condition (12.31) 

gives the velocity of this shock wave, which reads here: 



(3.2) At) = 4 (1 - *U *±« ,/ J- ^iStllM , with P(,) = i ,70 

2 \ q J 2 V q qo q-qo 70 
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Figure 1. The profiles depend on the friction coefficient 



We have the following result: 

Proposition 3.1. The shock wave propagates slower than the stain wave, and faster than the local 
wave speed c{qo{z)), that is 

Co qoizit))^ < z\t) < A . 



Proof: We fix t and set qs{t) and go = %{z{t))i as the left and right values of the shock wave. 
Then the shock velocity reads 



z'{t) = — ~ ^° ^ + ^0 / 1 Piqs) - Piqo) 



2 2 y q^qo g, - go 

We shall use some ^ &]qo,qs[ such that 



Pjqs) - Pjqo) 

qs - qo 



2cn 



P'iO = 



We have 



2 V ^sqo qs-qo V 2 2g^ 



where ^ = (i + ^j) • We get 



m < A ^gjU^) J^cge^ < A ■ 



2 2gJ V qsqo ' V2 2g. 



which reduces to 



qs , . ^ qref 
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and is equivalent to 

V W Qoizf) V go ' 

which is true since g^e/ > Qoi^f) and g^e/ > > ^ > Qo ■ 
On the other hand, since the expression 

A qs - go ^ Qs + go / 1 Pjqs) - Pjqo) 
2 gs 2 V ^'^^o ~ ^0 

is an increasing function of g^ for g^ > go, we get obviously the other inequality. (End of proof) 

Now we can construct the whole wave, made of a regular part corresponding to a part of the 
strain wave and a front shock whose position is determined by the Rankine Hugoniot condition (13.21) 
intertpreted as a differential equation whose solution z(t) gives the position of the shock. Figure 2 
shows the different positions of the shock wave for different values of the friction coefficient. 




Figure 2. The shock reduces the wave amplitude 



The velocity A is equal to 2932.5 m/s that is A = 1.78 Cq. The shock on the bottom of the 
ocean, at the depth zj = —3700m, corresponds to 15 times the value of the usual geostrophic 
pressure, thai is 

Pbottom = 363.6 10^ pascal = 363.6 bars , pref = 5.454 10*^ pascal = 5.454 kbars . 

This corresponds to strain densities 

qfond = 1.01288 , qref = 1.1296 . 

We notice that the variation from qbottom = Qoi^f) to g^e/ corresponds to an increasing of 11.5 % 
which is an increasing of about 1400 % of the pressure. The velocity of the wave is computed from 
the values 

Cref = Co q^°f = 2629.5 m/s , Wref = Cref —— = 303 m/s . 
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We notice that the velocities c^e/ and A = Wref + c^e/ are for more important than cq. The real 
profiles are drawn in red, and are shaped as a part of the strain wave profile cut by the slower shock 
wave. The eliminated parts are drawn in blue. We observe an important difference between the 
velocites of the different fronts, for several values of the friction coefficient. 

4. A WAVE OR NOT ? 

We denote by H (= —Zf here) the mean depth of the ocean. The wavelength A of the Saint- Venant 
waves must satisfy a condition of the form 

(4.1) A > 2 iV , 

where N the number of sonic interactions (back-and-forth) between the bottom and the surface of 
the ocean. This condition means that along a horizontal distance of a wavelength A, there are at 
least N such sonic interactions . The use of the Saint- Venant model is as more appropriate as 
is great. As in [2] for Rogue waves, we propose to require > 25, which implies for example a 
wavelength greater than 21400 m for an ocean depth of 3700m. A tsunami wave is expected when 
the condition (14.11) is fullfilled. A linear model for the surface elevation propagation was proposed 
in the historical paper [1] by K.Kajiura. 




Figure 3. The effect of the incidence 



The numerical tests in section 3 show the effect of the friction, whose effect is to shape up 
the amplitude of the wave near the front shock. At the same time, this front shock magnitude is 
eroded progressively as the wave propagates upwards to the surface. Since the amplitude of the 
front shock decreases, the velocity of this front shock decreases too. By following the propagation 
of a wave with an incidence of angle (that is only changing g into g coscj) and z into z / coscf)) , we 
get a longer path to travel with a weaker gravity constant) and a delayed wave compare to the case 
withou incidence. This is show on Figure 3. 

The value of the initial amplitude has a capital effect. It must be large enough to get, after 
erosion by the friction, a remaining wave near the surface which is sufficient to raise up the sea 
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surface and provoke a wave. The physical wave starts as a sperical wave, and propagates according 
to the incidences. The part with a small incidence will reach the surface later and will help the 
formation of the water wave. The part with a larger incidence will disappear because of the friction 
effect. The question of the value of the friction coefficient stays open, since the use of the strain 
density was never done before. It seems from the numerical tests that the correct values lay between 
1 and 10. Too large values provoke a sharp front wave which erodes rapidly and will never reach 
the surface with a sufficient amplitude to make a wave, which is not expected, since sometimes, 
tsunamis really occur. 

5. Annex: the source wave linearity 
We consider a general 2x2 hyperbolic system whose first equation has the form 

(5.1) qt + rrix = . 

We denote by Ai and A2 > Ai the eigenvalues of the flux matrix, which depend on q and m only. 
Then the general form of the second equation is 

m + (A1 + A2) - X1X2 Qx = S{q,m) , 

or 

(5.2) rrit + 2 u rUx + (c^ — m^) Qx = S{q,m) , 
by using the notations 

Ai + A2 A2 — Ai 

u = , c = . 

2 2 

and S{q, m) is a source term, assumed to be not identically zero. We have the following result: 

Theorem 5.1. The nonlinear non homogeneous system Ii5. li5.S\) admits non constant local 
solutions which are also solutions to the linear homogeneous system 

(5.3) qt + A q^ = , rut + A = , 

where A is a real constant, with the linkage Aq — m = B, another real constant. 

Proof: We look for local solution with a linkage of the form m = m{q). Then the system becomes 

qt + m'{q)q^ = , m{q) {qt + 2uq^) + [c^ - u^) qx = S{q,m{q)) . 

Since qt = —m\q) q^ , and S{q,m{q)) ^ 0, the second equation becomes 

{u{q, m{q)) - m'{q)f - c(g, m{q))^ _ 
S{q,m{q)) 

which has the form 

ij'{q) qx = I , 

{u{q, m(q)) — m\q))^ — c(g, m(q))^ 



by introducing a real function iplq) whose derivative is ip'lq) 

S{q,m{q)) 

Now, integrating with respect to x gives iplq) = x — K{t) , where K{t) is an integrating constant 
which may depend on t. Next, derivating with respect to t gives ip'iQ) Qt = — K'{t) , where we 
set qt = —m'{q) q^ .Hence we get m'(g) ip'{q) qx = K'(t) , and recalling that ip'{q) qx = 1 ,it 
remains m'{q) = K'{t). A new derivation with respect to x leads to m"{q) q^ = , and since the 
solution is not constant, we get m"{q) = 0, that is m{q) = Aq — B, with some constants A and B. 
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Next, K'{t) = m'{q) = A, and we have got tp{q) = x — At — xq , for some constant xq, or, locally, 

q = {x — Xq — At), which satisfies to (15.31) . since rrit + Anix = m'{q) {qt + Aqx) = .(end of 
proof) 

This result is a very general one, since no special hypotheses were needed on the second equa- 
tion (15.21) . Such waves are very common in the nature: water waves such as roll waves, rogue waves, 
tidal bore waves or also many other waves as reported in [2] or [4]. For example the double prop- 
erty of being either a solution to a non linear, non homogeneous systems and a linear homogeneous 
system provides the linkage between acoustics and gas dynamics in a wind instrument (see [3]). 

In the case of a conservative system, invariant by Galilean transform, the only choice of the 
function u is reduced to u{q, m) = m/q . We easily construct this way the usual Saint Venant system 
in hydraulics or the Euler equations in gas dynamics. 

In Section 3, the state at rest is not q = 0, m = 0, but q = qo, m = 0, so we look for a solution 
with the linkage At] — m = B since rj = q — qo = Q aX rest. We find a differentiel equation which is 
more complex than (15.41) . 
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